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I.  INTRODUCTION 


In  strategic  communication  systems  the  multiplicity  of  transmitters 
and  receivers  produces  a  large  ensemble  of  interference  sources  that 
must  be  carefully  examined  both  during  the  design  and  testing  phases.  To 
complicate  matters  further,  the  link  between  the  ith  transmitter  and  the 
jth  receiver  can  exhibit  nonlinear  characteristics,  thereby  generating 
intermodulation  and  cross-modulation  effects.  These  intermodulation  and 
crossmodulation  products  can  often  increase  to  intolerable  high  levels,  due  to 
accumulation  of  distortion  along  the  link.  Reduction,  then,  must  be  achieved 
through  compensation,  or  redesign  of  certain  subsystems.  Clearly,  in 
testing  such  communicatioii  systems  it  is  important  to  make  certain  that  the 
nonlinear  distortion  lies  within  acceptable  limits.  The  Volterra  series 
;  expansion  [l]-[3]  permits  description  of  a  nonlinear  system  in  a  compact 

form  and,  in  turn,  enables  computation  of  the  distortion  in  terms  of  the 
,  multivariable  transfer  functions  [A]-[6]. 

To  date,  however,  reliable  and  rapid  methods  for  finding  the  Volterra 
I  transfer  functions  (or  kernels)  from  laboratory  tests  have  been  lacking.  The 

work  by  Schetzen  requires  the  use  of  wideband  random  excitation  [7], 
which  is  somewhat  impractical  in  a  laboratory  environment,  both  because  suit- 
I  able  wideband  sources  with  guaranteed  uniform  spectral  density  in  the  region 

of  interest  are  difficult  to  find  and,  because  their  digital  analysis 
requires  unduly  high  sampling  rates.  Further,  the  method  requires  sufficient 
statistical  averaging,  involving  guesswork  on  the  part  of  the  test  engineer. 

’  Tlie  method  by  Weiner  and  Ewen  [  8]  uses  deterministic  input,  but  they  assume 

a  rather  restrictive  model,  namely  with  linear  transfer  function  H^(s)  and 
quadratic  transfer  function  H2(s^,S2)=3H^(s^)H^(s2)H^(Sj^+S2)  ,  where  'S'  is  a 
gain  constant.  In  realistic  systems,  the  modes  and  frequencies,  as  well 
as  the  residues,  of  the  quadratic  subsvstem  mav  be  different  from  those 

I 

suggested  by  this  model.  Further,  the  test  input  they  use  is  an  exponentially 
decaying  signal,  not  readily  available  from  standard  signal  sources. 

I  In  this  report  we  describe  a  method  for  determining  the  linear  and 

!  quadratic  subsystems  from  square-pulse  tests.  The  quadratic  subsystem  is 


1 


assumed  to  be  symmetric,  but  of  the  form  H-(8, ,s-)*H  (s,)H  (s,)H  (s,+s,) 
so  that  it  permits  sufficient  generality.  The  identification  method  involves 
two  transient  test  in  the  laboratory,  followed  by  analysis  by  the  computer. 

The  latter  consists  of  (a)  pole  determination  using  the  pencll-of-functlons 
method  [9]- [12]  and  (b)  computation  of  the  residues  by  a  least-squares 
technique.  Advantages  of  the  method  include  the  rapidity  of  the  laboratory 
tests,  as  contrasted  with  traditional  frequency-scan  approaches,  and  the 
explicit  determination  of  the  transfer  functions.  Furthermore,  the  method 
is  readily  extendible  to  H2(Sj^,S2,S2)  and  even  to  higher  order  transfer 
functions,  although  the  computations  grow  very  rapidly  for  these  cases. 

The  identification  technique  developed  here  represents  a  significant  improve¬ 
ment  over  existing  identification  techniques,  and  is  potentially  a  way  to 
turn-key  automatic  test  systems. 

The  structure  of  the  report  is  as  follows.  In  Section  II  we  study  the 
z-domaln  representation  of  quadratic  Volterra  subsystem.  Section  III 
discusses  a  general  method  for  computation  of  the  response  of  the  linear 
and  the  quadratic  subsystems.  In  Section  IV  the  bilinear  and  Volterra  responses 
of  the  quadratic  subsystem  to  square  pulses  is  computed  in  the  form  of 
explici:  formulas.  The  form  of  the  responses  shows  they  are  representable 
as  the  step,  or  impulse  responses  of  certain  equivalent  linear  systems. 

As  a  consequence,  the  pencil-of-functions  method  by  Jain [10]- [12]  can 
be  used  to  determine  the  poles  of  the  quadratic  subsystem.  The  complete 
identification  procedure  for  a  two-variable  system  is  discussed  in  Section 
V.  Computer  generated  examples,  also  presented  in  Section  V,  demonstrate 
the  success  of  the  approach  developed  in  this  report. 
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II.  z-DOMAIN  CHARACTERIZATION  OF  THE  QUADRATIC  VOLTERRA  SYSTEM 

In  this  section  we  shall  be  Interested  in  the  z-domaln  characterization 
of  the  quadratic  nonlinear  system  shown  in  Fig.  1.  This  system  is  the  most 
basic  quadratic  Volterra  system,  henceforth  referred  to  as  the  single  multi¬ 
plier  (IM)  configuration.  It  will  be  shown  that  this  system  can  be  repre¬ 
sented  mathematically  by  the  second  degree  term  in  the  Volterra  series  [14], 

[  7]  with  the  corresponding  quadratic  transfer  function  H  (s  ,  s„)  *  H  (s  ) 

1.  ^  SL  A. 

H  (s„)H  (s,+  S-).  For  the  case  where  H  ,  H,  and  H  are  rational,  we  will 
b2cl2  a  bc 

find  the  z-domaln  description  of  this  system  for  impulse  and  step  Invariant 
criteria.  Throughout  we  will  assume  these  components  are  causal  and  linear. 


x(t) 


y2(t) 


Fig.  1.  Basic  quadratic  Volterra  system 

2. 1  Continuous-Time  Analysis 

From  Fig.  1  we  observe  the  output  of  the  multiplier  to  be 


x^(t) 


h  (C)  x(t-5)  dK 
a 


h^(ri)  x(t-ri)  dn 


(1) 


From  this  the  response  of  block  H^,  i.e.,  the  system  output,  is  found  to  be 


y2(t) 


h  (t)  x  (t-T)  dx 
c  c 


h  (X) 
c 


h  (O  xit-T-O  dK 
a 


h^(n)  x(t-x-ri)  dn  dx 


(2) 
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Letting  -  T  +  and  T2  *  T  +  n,  we  have 


y2(t) 


^^’*^2'’’^^  x(t-T^)  x(t-T2)  dT^  dT2  dT 


r  r 


x(t-Ti)  x(t-T2)  dT2 


(3) 


where 


h^(T)  h^(T^-T)  h^(T2-T)  dX 


(A) 


The  two-dimensional  Laplace-transfonn  of  h2  is 


-(s,T^+  s,T_) 

hc(T)  h^(T2-T)  e  dXj^  dT2  dx 


■(Si+  So)X 


5iX-  s,/r  r  -s,  (X,-  x) 

h^(T)  e  ^  2  J  h^(x^-  X)  e  ^  ^  dx^ 


hj^(X2-  X)  e 


dX2  dx 


=  H^(8^)H^^(s2)H^(s^+  S2) 

We  will  omit  the  qualifier  "two-dimenSlonal"  when  it  is  clear  from  the  context. 
However,  the  significance  of  this  qualifier  becomes  particularly  evident  when 


we  associate  with  y2(t)  in  (3)  the  following  two-dimensional  response 


y2(tj^,t2)  =  h2(Tj^,X2)  x(tj^-  x^^)  x(t2-  T2)  dx^^  dX2 


(6) 


Straightforward  application  of  the  definition  of  two-dimensional  Laplace 
transform  and  some  manipulation,  yields 
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S2)  X(s^)  X{S2) 


(7a) 


S2) 


«  H2(Sj,  S2)  X(Sj,  S2) 


(7b) 


which  represents  the  input-output  relationship  of  a  two-dimensional  linear 
system  as  shown  in  Fig.  2.  We  will  call  associated  two- 

dimensional  response  since  y2(t)  *  y2^*'l’  *'2^1 


x(tj^,t2) 


Fig.  2.  Associated  two-dimensional  linear  system 


Rational  Case 

When 

H  - 
a 

n 

1 

i=l 

A. 

1 

s+a. 

1 

n 

B. 

s+b . 

J 

»b  - 

Z 

j-1 

n 

H 

c 

1 

£-1 

or,  equivalently 

9 

n 

h^(t) 

Z  A. 
i=l  ^ 

n 

h^(t) 

* 

Z  B. 

n  -c  t 

h  (t)  =  ^  C  ^ 

d-l 


t  >  0 


(8a) 


(8b) 
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then  equations  (3),  (6)  and  (7)  may  be  written  as 


00 

,(t)  -  z  I  "'■2^  ‘^’'^1  ‘^’^2 

_oo 


y^(t^,t^)  «  I  JJ  x(t^-  T^)  xCt^-  T^)  dT^  dT^  (9b) 

1  «  I  •  X»*  X 

*  '  rr\ 


Y^(s^,s^)  -  I  ^ij£  X(s^,S2) 

1 1  j  j  1 


where 


(s^+a^)(s2+b^)(sj^+s2+cp 


(10a) 


(10b) 


It  follows  from  (9)- (10)  that  the  response  is  a  weighted  triple  sum  over 
the  response  of  the  elementary  quadratic  system  shown  in  Fig.  3  where  a  =  a^, 
b  =  bj ,  and  c  »  c^.  We  will  abbreviate  this  process  of  triple  summation  as  WTS. 
Thus,  we  arrive  at  the  important  conclusion  that  it  is  only  necessary  to  analyze 
the  elementary  system  of  Fig.  3  in  detail;  the  properties  of  the  original  system, 
e.g.,  its  response,  would  follow  immediately  through  the  appropriate  WTS  process. 


y2^t) 


Fig.  3.  Elementary  quadratic  system 


Impulse  response  of  the  elementary  quadratic  system 


Using  h  (t)  =  e  u(t) ,  a  =  a,b,c  in  equation  (4)  we  obtain  the  quadratic 


J 


Volterra  kernel  for  the  elementary  system  of  Fig.  3  as 


gCXj^.T^)  -  I 


_CT  -a(T^-T)  -b(T2-  T) 

e  u(t)  e  u(t^-  t)  e  uCt^-  t)  dx 


-(c-a-b)T  _  .  "®^1  "^^2 

e  dx  u(t*)  e  e 


-ax  -bx.  -(c-a-b)x* 

=  e  e  ^  ^  -  u(T.)  u(T-)  (U) 

c— a-b  1  z 

where  T*  =  MinCXj^,  ^2) ,  it  is  assumed  that  c  a  +  b. 

The  impulse  response  of  the  elementary  quadratic  system  is  obtained  by 
setting  “  T2  =  t  in  (11): 


y2(t) 


Example  1 


-(a+b)t  -ct 
e  —  e 

c  -  a  -  b 


For  a  =  1,  b  =  3  and  c  »  6  we  have 


-T.  -  3T-  -  _  -2t* 

g(Ti,T2)  =  e  - 1 -  u(Tj^)  u(T2) 


This  is  depicted  In  Fig.  A  where  the  hatched  region  denotes  the  region  of 
zero  value  for  g(Tj^,T2). 


Note  also  that  the  two-dimensional  response  to  the  input 


x(t^,t2)  “  <5(t^-a)  6(t2-S) 


y2(tj^,t2)  =  e 


-a(t^-a)  -b(t2-  8)  ^  _  ^-(c  -  a  -  b)t* 


c  -  a  -  b 


u(t^  -  a)  u(t2  -  6) 


The  upper  limit  x*  implies  that  integration  is  done  from  0  to  or  X2. 
ever  is  cnnller. 


whieh- 
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typical  quadratic  Kernel 


where  t*  «  min  -  a,  t2  -  Bl.  This  follows  from  (11)  and  the"8tationarit3)' 
of  the  associated  two  dimensional  system. 


Bilinear  Operator  [7 ] 

Equation  (3)  may  be  rewritten  as 

y2(t)  =  V2[x(t)] 


J  h2(Tj^,T2)  x(t-Tj^)  x(t-T2)  dTj  dT2 


That  is  y2(t)  may  be  viewed  as  the  result  of  a  nonlinear  operator  upon  the 

input  signal  x(t).  This  operator  may  be  generalized  as 

00 

V2{xj^(t),,X2(t)}  *  II  h2(Tj^,T2)  dr^  (13) 


Clearly,  V2{->*}  is  bilinear  in  and  X2;  we  will  call  it  the  bilinear 
Volterra  operator,  or  simply  the  bilinear  operator.  Note  that 
V2  [  X  ]  =  V2{x,  x} 

The  significance  of  the  bilinear  operator  arises  from  the  fact  that 

n  n 
j 

Example  2(a) 

We  will  find  the  bilinear  response  of  the  elementary  system  of  Fig.  3 
for  the  case 

x^(t)  =  6(t-a),  X2(t)  »  6(t-8) 


V  [  Z  a  X  J 
i=l 


Kl  It 

E  E  a  a,  V„{x  ,  x  } 

i-1  j-1  ^  J  J 


(14) 


where  we  will  assume  a  <  6. 


From  equation  (15)  we  find  the  bilinear  response  to  be  the  convolution 

of  h2(tj^,t2)  and  6(tj^-a,  t2-3) ,  with  tj^  and  t2  set  equal  to  t.  This  is 
equivalent  to  the  integral  of  the  product  of  h2('r2^»T2)  ^  two-dimensional 

impulse  located  at  =  t  -  a  and  T2  =  t  -  3.  This  is  depicted  in  Fig.  5a  for 
four  different  values  of  t,  namely  t  =  0,  a,  3  and  t  >  3.  Now,  for  the 
particular  case  of  the  1  -  M  system  of  Fig.  3,  we  have  (denoting  l/(c-a-b)  as  (J)) 


V2{xi. 


X2} 


00  J 

■r  1  -aT^  -bT2 
e  e 

00 


-(c-a-b)T^ 

1-e _ 

c  -  a  -  b 


u(T2)  6(t-3-T2)  dT2  u(Tj^)  6(t-a-T^)  dr^ 
-aT 


(5(t-a-T^)|  [e  "  -e 


-bT2  -(c-a)T2 


]  6(t-3-T2)dT2  dT^ 


r°°  -at 


^  6(t-a-T^)le  _g-(c-a)(t-3)^  u(T^-t+3)dT^ 


^^-b(t-3)  _^-(c-a)(t-3)j  g  ^"^1  6(t-a-T^)dT^  u(t-3) 


t-3 


»  (f.  _^-(c-a)(t-B)j 

.  1  ^-a(t-a)^^-b(t-3)  _^-(c-a)(t-3)^ 

c  -  a  -  b 

For  clarification  of  limits  of  integrations  leading  to  steps  two  and  three 
see  Fig. 5b.  A  similar  formula  can  be  derived  for  the  case  a  >  3. 
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tr> 


■  « 


(15) 


Note  that  the  transform  of  is 

V2{X^(s),  X^Cs)}  =  H2(Sj^,82)X^(Sj^)X2(82) 

and,  clearly,  (15)  can  also  be  written  in  the  transform  domain: 


V  [  S  a  X  (s)]  = 

i=l 


It 

Z  Z  a  a  V  {X  (s),  X  (s); 
1=1  j=l  ^  J  ^  J 


To  demonstrate  the  usefulness  of  the  transform  domain  we  consider  Example 
2(a)  again. 

Example  2(b) 


Since 


X^(Sj^)  =  e 


X2(S2)  =  e 


(s^+a)  (s2+h)  (Sj^+S2+c) 


-as^  -Bs2 
e  e 


Application  of  George's  theorem  (Corollary-Appendix  A''  gives  immediately  the  same 
time  domain  bilinear  response  as  before. 

Finally,  before  leaving  the  discussion  on  bilinear  response,  we  note  that 
the  bilinear  response  of  the  general  1-M  section  with  H^,  H^,  and  characterized 
by  (8),  is 

n 

V2{X^,X2}  =  Z  X  (s  )X  (s  )  (17) 

The  time-domain  response  is  of  course  the  inverse  transform  of  (17). 

Bilinear  Response  From  Block  Diagram  - 

We  will  show  that  the  bilinear  response  of  the  1-M  (Fig.  1)  section  can 
be  obtained  by  applying  Xj^(t)  to  H^(s)  and  X2(t)  to  H^(s)  as  shown  in  Fig.  6(b). 

Of  course,  this  procedure  is  useful  only  for  pencil-and-paper  purposes  and 
simulation;  it  is  not  useful  for  practical  situations  where  it  is  generally 
not  possible  to  Isolate  H^(s)  and  Hj^(s). 


(a)  Block-d iaqram  for  bilinear  response 


X,  (t) 


(b)  Symbolic  representation  of  bilinear  response 
Fig.  6.  Bilinear  response  from  block  diagram 
Substituting  for  h2(Tj^,T2)  from  (4)  and  interchanging  the  order  of 
integration,  we  have 


KJO 

hc(T) 

J 

h^(Tj^-T)x^(t-T^)dT^ 

. 

h^(T2-T)x2(t-T2)dT2dT 


X*!— 00 


or,  with  -  T  and  =  T2  -T, 


»00  -00  *00 

V2{xi,X2}  =  h^(T)  h^(?^)x^(t-T-C^)dS^  J  h^(C2)x2(t-T-52)dC2  dT 


js— 00  —00 


hc(T)  YgCt-T)  yj,(t-T)  dt 


,ro 


I 

J  c 


h  (t)  X  (t-T)  dr 


where  yg(t)  is  the  response  of  H^(s)  to  the  input  Xj^(t) ,  y^^Ct)  the  response 

of  li  (s)  to  the  input  x.,(t),  and  x  (t)  =  y  (t)y,  (t).  Therefore,  the  inter- 
D  z  c  a  D 

pretation  of  (15)  in  the  form  of  Figi  6(a),  and  Fig.. -6(b),  is  justified. 


The  concept  of  impulse  invariance  of  a  two-dimensional  linear  system  can 


be  enunciated  as  follows.  Given  that  a  continuous  input  signal  xCtj^,  12)  is 
applied  to  the  sampled-data  system  of  Fig.  7a;  then  the  system  of  Fig.  7b  is 
said  to  be  impulse-invariant  if  its  response  y(k^,  k.2)  to  the  sequence 
x(k^,  k2)  =  x(k^A,  k2A)  equals  the  sampled  response  y(k^A,  k2A)  of  Fig.  7a. 


x*(t^,t2) 


(tj^,t2) 


Ideal 

Sampler 


y(k^A,k2A) 


x(k^,k2) 


y(kj^,k2') 


Sampler 


x*(t^,t2)  *  SE  x(kj^A,k2A)6(tj^-kA^,t2-kA2) 

k  k  00 

*^1  *  *^2 

(a)  (b) 

Fig.  7.  Impulse  invariant  discrete-time  system 


Note  that  the  above  concept  of  impulse-invariance  is  a  direct  extension 
of  this  concept  for  the  one-dimensional  case  [13]. 


We  will  use  the  known  result  that  the  impulse  invariant  equivalent 
z-domain  function  is  given  by  the  formula  [14] 


FCVi.Vz)  dv^  dv2 

-(s.-vJA  -(s,-v„)A 
(1-e  ^  ^  )(l-e  ^  ^  ) 

Specifically,  for  the  function  G(Sj^,S2)  “  l/(s^+  a)(s2+  b)(sj^+  S2+  c)  , 
the  quadratic  transfer  function  of  Fig.  3,  we  have 


F(z^,Z2) 


+i< 


s .  A 
1 


^2TTj^ 


(18) 


«h«re 

Q(« 

Then, 

G(z 


G(Zj, 


z^-e 


SiA 


'2irj^ 


-J 


dVi  dvj 


-fe  -V, )A  -(s, -v^A 
(Vj^+a)(v2-H>)(v^+V2+c)(l-e  ^  ^  )(l-e  ^ 


■<sj>| 


dv. 


-J®  (v2+b)(l-e 

-(s  +a)A  j 

i-  ••  ( 


-(82-V2)A 


-(8.+a)A  -(a  +v,+c)A 

(v2+c-a)  1-e  1-e 


1-e 


_ L.  f 

-(Sj+a)!  2iTj  J 


1-e 


-(v2+c-a)A 


(82-V2)A 


-j"  (v2-H>)  (v2+c-a)  (1-e  )  (1-e 


■(81+V2+CJA  ‘^''2 


-(8.+a)A  j® 

e _  f  1-e 

-(8,+a)A  2nj  J 


-(v2+c-a)A 


dv- 


1-e  -J®  (v2+b)(v2+c-a)(l-e 


-(8,+v,+c)A  -(s.-v  )A 

)  (1-e  ^  ) 


-(Bj^+a)A 


e _  _  1  J  2irm, 

-(8.+a)A  ^  A  ^^®2  ^  A  ^ 

-  1  Bl*— 

1-e 


-(8j+a)A 

-(e,+a)A  -(8,+8,+c)A  _  I 


(1-e  ^  )(l-e  ) 


.)  - 


1-e 


P(e2) 

-Ta^+i^+cjA 


1-e 


-(82+c-a)A 


-(8^+8„+c)A 
(e2+b) (82+c-a) (1-e  ) 


(20) 


1,*2> 


I-«^T'-i77rT^cA7^T~3i, 


(^-bA  _^-(c-a)Aj  ^^-1 


(c-a-b)(l-e-*^  zf")(l-e-‘=^^-lz2'b  (l-e-("-“>^  Z2-^)  d-e"'’^ 


-1  -1  / 
*1  *2 


(c-b-a)(l-p  z^~  )(l_q  z^'^Xl-r  Z2"^) 

with  p  ■  e”*^  ,  q  -  ,  and  r  - 


(21) 
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Note  in  the  above  we  evaluate  line  integrals  by  contour  integration,  closing 
for  convenience,  in  the  left  half  plane  in  the  first  integral  and  in  the  right 


half  plane  in  the  second.  Also  note  Re[Sj^]  >  0,  Re[s2]  >  0  since  h2('rj^,T2)  is 
causal . 

example  3 

For  a  =  1.025866,  b  *  3.250379,  c  *  5.753641  and  A  »  0.05  s,  we  have 

^  0.0575z^~^ 

(s^+a)  (s2+b)(Sj^+S2+c)  l .  477397(1-0. 95z^"^)  (1-0. 8522"^)  (1-0.  Z2~^) 

Because  of  its  Importance  we  give  the  conversion  pair,  derived  above, 
explicitly : 


(pq-r) 


(s^+a) (s2+b) (s^+S2+c) 


(c-a-b)  d-pz^'^)  (l-qz2"^)  d-rzj”^  Z2 


(23) 


(note  p,  q  and  r  are  the  z-plane  maps  of  the  poles  -a,  ~b  and  -c;  e.g., 
-aA 

p  -  e  ,) 


An  alternative  derivation  (in  the  time  domain)  of  the  above  conversion 
pair  can  also  be  given.  A  block  diagram  description  of  the  input-output 

relationship  (of  the  elementary  quadratic  system  of  Fig.  3)  in  the  z-domain 
is  given  in  Fig.  8. 


x(k) 


y(k) 


a=  -e~‘^^)/(c-a-b) 

Fig.  8.  Impulse-invariant  z-domain  characterization  of 
the  elementary  quadratic  system 
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Examp le  4 


■k 


Let  us  calculate  the  response  of  the  systenn  of  Example  2  for  four  sample 
points  with  k  “  0,  1,  2,  3  and  the  input  1,  1,  0,  0,  ...  We  find 

y  (k):  0  1  1.95  1.8525 

a 

y^(k):  0  1  1.85  1.5725 

so  that  their  product  gives  the  following  input  to  the  output  block 

X  (k):  0  1  3.6075  2.91306 

c 

Then  y(k),  as  obtained  from  the  difference  equation  y(k)  “  0.75y(k-l)  + 

01^  X  (k)  ,  a  *  0.03892,  is 
c 

y(k):  0  0.00151  0.00660  0.00936 

STEP- INVARIANT  z-TRANSFORM  FOR  1-M  SECTION 

The  concept  of  step  Invariance  of  a  two-dimensional  linear  system  can  be 
enunciated  as  follows.  Given  that  a  continuous  input  signal  x(t^,  t^)  is 
applied  to  the  sampled-data  system  of  Fig.  9a;  then  the  system  of  Fig.  9b  is 
said  to  be  step-invariant  if  its  response  y(k^,  k2)  to  the  sequence  x(k^,  k2)  = 
xCk^A.k^A)  equals  the  sampled  response  y(k^A,k2A.)  ■  of  Fig.  9a. 


p(t)  »  1  for  0  _<  t  _<  A,  0  otherwise 

(a)  sampled-data  system 
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x(k^,k2) 


y(kj^,k2) 


(b)  Digital  system 

Fig.  9.  Step- invariant  digital  system 

The  above  concept  immediately  results  in  the  following  definition. 

Definition.  Given  a  continuous  two-dimensional  system  H  .  a  dieital  twn- 

a 

dimensional  system  is  said  to  be  step  invariant  if  the  response  of  to 
u(k^)  u(k2)  equals  the  response  of  to  u(t^)  u(t2)  at  the  sampled  points. 

Note  that  the  above  concept  of  step-invariance  is  a  direct  extension 
of  this  concept  for  the  one  dimensional  case  [13]. 

We  begin  the  derivation  by  observing  that  the  transfer  function  of  the 
zero  order  hold  is 


\old^®l’ 


(1 


-s^A  -s„A 

e  M(1  -  e  ^  ) 


(24) 


and  the  transfer  function  of  the  system  between  the  ideal  samplers  (in  Fig. 
9a)  is 


% 


-s.  A 

(1  -  e  ^  )(1 


-s-A 
e  ) 


®1®2 


-s^A  -S2A  -(sj  +  S2)A 
(1  -  e _ -  e _ +  e _ )_ 

"l"2 


H2(s^, 


S2) 


-  H2(s^,  S2)  + 


(25.1) 


_ 1_ 

®1®2 
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-s^A 

- -  H  (s  ,  s_) 

Sj^S2  2  12 


-S2A 


s^S2 


-(Sj^  +  s^)h 


"l"2 


^2^®1’  ®2^ 


(25.2) 


(25.3) 


(25.4) 


Further,  we  make  the  Important  observation  that  the  step-invariance  of 
H2(s^,  S2)  to  some  H2(Zj^,  z^)  is  equivalent  to  the  Impulse  invariance  of 
^^2(3^,  S2)  to  the  aameH2(z2^»  .  This  is  obvious  from  the  enunciation  of  the 
basic  concepts  of  impulse  and  step  Invariance. 


We  will  now  restrict  our  attention  to  the  IM  section  so  that 

„  ,  .  ABC  _ _ 

2^®1’  ^2"^^  (Si+  a)  (s»+  b)(sT+  s-+  c) 


(26) 


■'2^'  *^2 

Observe  that  each  of  the  four  terms  on  the  right  hand  side  of  (25)  is  of  the  form 


•(ysj^  +  VS2) 


^2,i  ^®1’®2^ 


ABC 


"l"2 


(s^+  a)(s2+  b)(sj^+  S2+  c) 


y,  V  •=  0  or  A 

i  -  1,  2,  3,  4 

Partial  fraction  expansion  yields 


(27) 


,  .ABC  -(Vis^+  VS2)  1 


_  _i_]  _ I _ 

s^+a  S2  S2+b  (Sj^+  S2+  c) 


ABC  1 

— r  e  [ - 

ab  s,s 


j^„2  s^(s2+  b)  ^®1^  a)(s2+  b) 


Si+  S2+  c 


(28) 
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In  turn,  each  of  the  four  terms  on  the  right  hand  side  of  (29)  is  of  the  form 

(29) 


G(s^,S2)  •  K  e 


■(ViSj^+  \)s^) 


(s^+  a)(s2+  B)(Sj+  s^+y) 


for  which  the  impulse-invariant  z-transform  (from  (23))  is 

-1  -1, 


G(z^,22)  -  K  Zj  ^  Z2  ^ 


Zf  (iT^-p) 


(Y-a-B)  (1-TtZj^  ^)(1-^Z2  ^)(l-pz^  ^Z2  ^) 


(30) 


Using  the  property  of  linearity,  the  step-invariant  z-transform  of 
H2(s^,S2)  consists  of  sixteen  terms  obtained  by  using  the  partial 
fraction  expansion  of  (28)  in  (25)  and  finally  using  the  correspondence 
of  (29)- (30)  upon  each  term.  The  result  is 


^2^^1’^2^ 


-1  -1 

ABC  ^1  ^2 

ab 


-1  -1 
r  z^  Z2 


< 


1  -  r  +  I 


1  -  z 


-1 


c  -  a  ,  -1 

V  1  -  P 


1  -  z~^ 

+  3ZL  - 2 - 

c-b  ,  -1 

1  -  q  z. 


,  -aA 

where  p  =  e  ,  q 


-bA  -cA 

e  ,  r  “  e 


+  P<l-r 
c-a-b 


(l-z^"^)(l-Z2"^) 
(1-pz^  ^)(l-q  Z2  ^) 


This  is  represented  in  block-diagram  form  in  Fig. 10. 
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SQUARE R 


Block  diagram  of  s tep- i nvar i ant  z-transform 
for  the  1  -  M  sect i on 


III.  DFT  BASED  COMPUTATION  OF  QUADRATIC  VOLTERRA  RESPONSE 


A  general  method  for  computation  of  the  response  of  a  second  order 
Volterra  system  (see  Fig.  11)  is  highly  useful  in  the  investigation  of 
Volterra  system  identification.  The  features  desired  in  such  a  method  are 
the  following: 

a.  The  characterization  of  the  blocks  in  Fig.  11  should  be  permitted 
in 


s  domain 
z  domain 

f  domain  (i.e.,  in  terms  of  the  frequency  ch.  H(f(k))) 
t  domain  (i.e.,  in  terms  of  the  impulse  response  h(t(k))) 

b.  The  specification  of  the  input  should  be  permissible  in  the  time 
domain.  This  is  useful  in  transient  type  of  tests,  such  as  with 
pulse  Inputs. 

c.  The  specification  of  the  input  should  be  permissible  in  the  frequency 
domain.  This  is  useful  in  steady-state  sinusoidal  tests. 


Fig.  11.  Quadratic  Volterra  system 
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Such  a  method  can  be  best  developed  by  use  of  the  discrete  Fourier 
transform  (DFT) .  We  shall  In  fai-t  use  the  fast,  algorithmic  version  commonly 
known  as  the  fast  Fourier  transform  (FFT) .  And,  we  will  use  the  terms  DFT 
and  FFT  interchangably .  To  develop  the  method  we  begin  first  with  the  case 
of  the  linear  system. 


3.1  Linear  System  Response  via  DFT 

Consider  a  linear  system  described  by  a  z-domaln  transfer  function  H(z). 

Note  that  if  an  s-domain  transfer  function  is  given  instead,  then  its  z-domain 
pulse-equivalent  function  can  be  obtained  by  use  of  the  FORTRAN  program  STOZllO], 
An  example  is  the  following 


(lO^)s 

s^  +  4(10^)s  +  10® 


0.02z~^  -  0.0196117Z  ^ 

1  -  1.8848022"^  +  0.923116z”^ 


where  a  sampling  of  A  =  0.02  ms.  has  been  used. 

The  choice  of  the  sampling  interval  and  the  number  of  points  to  be  used 
in  simulation  generally  requires  some  experimentation.  However,  certain 
guidelines  may  be  given  based  upon  signal  theoretic  considerations.  These 
will  be  discussed  later,  and  we  assume  that  the  sampling  interval  A,  an  input 
(of  essential  length  b  second,  or  B  samples)  x(t),  and  the  number  of  points 
N  to  be  used  for  DFT  processing  have  been  selected. 


Notation 


A  =  Sampling  interval 
N  =  Number  of  samples  used  in  DFT 


L  =  ^  +  1 


T  =  NA 


6  =  i.e.,  spacing  of  DFT  frequencies, 

or  frequency  resolution  in  DFT 
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H(s)  =  s-doraaln  transfer  function 


H(z)  =  z-domain  transfer  function 

H(f^)  =  Frequency  characteristic  of  the  network, 

where  f  =  m 
m 

H(tj^)  =  Impulse  response  of  the  network,  where  tj^  "  kA 


f  =  L6  =  (-^  +  l)iS,  i.e.,  the  highest  unambiguous  frequency 

II13X  z. 

Whatever  the  specification  of  the  linear  system,  we  convert  it  to  the 

H(f  )  form  as  shown  in  the  flowchart  of  Fig.  12. 
m 


m 


Fig.  12.  Conversion  of  system  specification  to  conjugate 
symmetric  H(f  )  form 


Also,  the  input  is  converted  to  the  frequency  form  as  shown  in  Fig.  13. 


Time  pulse  x(t,  )  X(f  ) 

K  m 


X(f  )  form, 
m 


For  purposes  of  quick  reference  we  give  below  the  pair  of  discrete  Fourier 
transform  formulas  : 

N-1  -j  mk 

X(f  )  =  Z  x(t,  )  e  ^  (32a) 

k=0  ^ 


x(t^) 


1 

N 


N-1 

E 

m=0 


X(f  )  e 
m 


2it 


(32b) 


These  relationships  will  also  be  written  for  convenience  as 


x(tk) 


DFT  . 
> 

dft“ 


X(f  ) 
m 
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Then  the  network  response  can  be  obtained  as 


Y(f  )  =  H(f  )  X(f  ) 
m  mm 

y(t.  )  =  DFT“^  Y(f  ) 

K  m 

We  remark  for  completeness  that  y(t^), 
the  following  circular  convolution  [3] 


(33a) 

(33b) 

the  result  of  (33b)  is  equivalent  to 


y(t,) 


N-l 

I  x(t  )  h((t  -  t  )) 
j=0  3  *^3 


(34) 


where  ((t,  -  t.))  denotes  the  difference  of  t,  and  t .  on  a  circular  (or 
k  1  ‘^1 

periodic)  basis  on  a  circle  of  circumference  T  =  NA.  Formula  (34)  brings  into 
focus  two  important  facts.  Formula  (34),  or  equivalently  (33b),  can  approxi¬ 
mate  the  familiar  aperiodic  convolution  integral 


y(t) 


f  <50 

x(t)  h(t-T)  dx 

«  00 


if 


a.  it  is  multiplied  by  A,  and 

b.  A+B-lj<N,  ora  +  b<T 
where 

a  =  duration  of  impulse  response  in  seconds 

A  =  duration  of  impulse  response  in  samples 

b  =  duration  of  input  pulse  in  seconds 

B  =  duration  of  input  pulse  in  samples 


(35) 


(36) 
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Example  5 

The  simple  first  order  system  considered  is 

1  ^  0.09754z~^ 

®  1  -  0.95123z"^ 

which  has  a  cutoff  frequency  at  approximately  0.08  Hz.  Since  the  roll-off 
rate  is  -20dB/decade,  it  is  reasonable  to  take  the  highest  frequency  of  in¬ 
terest  as  f  =  5  Hz.  Thus  we  take  the  sampling  interval  to  be  A  =  0.1s. 

Now  since  the  time  constant  of  the  network  is  2  sec,  we  should  estimate  the 
Impulse  response  to  be  about  8s.  Allowing  for  an  input  duration  of  Is,  we 
could  take  T  to  be  10s.  However,  for  reasons  of  frequency  resolution  (dis¬ 
cussed  later)  we  use 
N  =  256 

T  =  NA  =  25.6s 

Two  test  inputs  are  examined  below. 

(a)  Square  pulse  input  -  width  b  =  1  sec  (B  =  10  samples) 

The  true  response  is  easily  found  to  be 

y(t)  =  2[1  -  u(t)  -  2  [1  -  u(t  -  1) 

where  u(t)  denotes  the  unit  step  function.  The  samples  y(tj^)  computed  by  the 
DFT  method  (via  program  VOL  2  -  Section  IV)  are  compared  with  the  true  value 
in  Table  1;  a  graphical  comparison  is  given  in  Fig.  I4. 


Table  1 


0 

0.1 

1.0 

2.0 

5.0 

True  y(tj^) 

0.0 

0.442 

0.107 

DFT  based  y(tj^) 

0.0 

0.0975 

0.190 

0.477 

0.106 
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RMPLITUDE 


Fig.  14.  Comparison  of  true  and  DFT-based  responses 
(Example  5  -  Square  pulse  input) 


Exponential  Input  -  x(t)  =  e 

Because  of  the  long  duration  of  the  input  the  number  of  points  was  in¬ 
creased  to  N  =  512. 

The  true  response  is  readily  checked  to  be 


,,-0.25t  -O.St.  f  ^ 
y(t)  =  4(e  -  e  )  u(t) 


Tlie  samples  y(t,  )  computed  by  the  DFT  method  (via  program  V0L2)  are  compared 
K 

with  the  true  values  in  Table  2;  a  graphical  comparison  is  given  in  Fig.  15. 


Table  2 


0.2 

2.0 

5.0 

10,0 

True  y(tj^) 

0.415 

- 1 

0.689 

0.818 

0.301 

DFT  based  y(t,  ) 
k 

0.0 

0.0975 

0. 188 

0.420 

0.698 

0.967 

0.828 

0.305 

T* 


C  JO  C  (>C  0  JC  I.  I?0  jot  it!  211  2‘1  JK  )00 


1 TME 

Fig.  15.  Comparison  of  true  and  DFT-based  responses 
(Example  5  ~  Exponential  input) 

It  is  seen  above  that  the  DFT  method  is  able  to  compute  the  network 
response  quite  accurately.  For  further  accuracy,  the  sampling  interval  must 
be  decreased  or  the  number  of  points  must  be  increased  or  both.  We  should 
also  remark  that  if  the  input  were  e  a  sampling  interval  of  0.05s.  would 
be  used  rather  than  0.1s. 

Example  ^ 

The  oscillatory  system  considered  is 

_ (10‘^)3 _  ^  0.02z~^  -  0.0196117z~^ 

s^  +  4(10^)s  +10®  1  -  1.884802z"^  +  0.923116z“^ 

which  has  a  natural  frequency  at  approximately  1.6  KHz.  Since  the  roll-off 

rate  is  -40dB/decade  it  is  reasonable  to  take  the  highest  frequency  of  interest 


as  f  =25  KHz.  Thus  we  take  the  sampling  interval  to  be  A  =  0.02  ms. 
max 
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Further,  since  the  time  constant  is  0.5  ms,  we  should  estimate  the  duration  of 
the  impulse  response  to  be  about  2  ms.  Allowing  for  an  input  signal  duration 
of  0.5  ms,  we  could  take  T  to  be  2.5  ms.  However,  for  reasons  of  frequency 
resolution  (discussed  later)  we  use 

N  =  512 

T  =  NA  =  10.24  ms 

The  test  input  is  taken  to  be  a  square  pulse  of  width  b  =  0.6  ms.  (B  = 

30  samples).  Correspondingly , the  true  response  is 

y(t)  =  sin(lO^t)  u(t)  -  ^  sln(10^ (t-b) )  u(t-b). 

The  samples  y(t^)  computed  by  the  DFT  method  (via  program  V0L2)  are  compared 
with  the  tru  values  in  Table  3;  a  graphical  comparison  is  given  in  Fig.  16  , 


Table  3 


t ,  ms 
k 

0 

0.02 

0.04  1 

1 

0.1 

0.2 

0.5 

1.0 

True  y(tj^) 

0.359  ’ 

0.689 

1 

0.610 

-0.353 

0.266 

DFT  based  y(tj^) 

0.0 

0.191 

0.360 

0.694 

0.633 

-0.369 

0.271 

In  this  final  example  the  network  specification  is  given  in  terms  of  the 
frequency  response 

H(f)  =  1  for  |f|  <  5  KHz 

and  a  square  pulse  input 

x(t)  =  1  for  0  ^  t  0.2  ms 

Since  the  input  is  a  square  pulse  of  width  0.2  ms  the  highest  frequency 

of  interest  f  is  estimated  at  10  f  =  50  KHz.  Thus  the  sampling  frequency 

max  zero  o  m  / 

must  be  lOO  KHz  requiring  A  =10iis.  Then  we  have 
Input  duration  =  0.2  ms  (20  samples) 

Impulse  response 

duration  =  20  t  =  2  ms 

zero 

=  200  samples 

We  therefore  take  N  =  512.  The  true  response  is  known  to  be 
y(t)  =  ^  {Si[107T(10^)t]  -  Si[10Tr(10^)(t  -  0.2(10“^))]} 

where  ^ 

Si[v]  =  j  li^dv 

o 

The  samples  computed  by  the  DFT  method  are  compared  with  the  true  value 

in  Table  3.  A  graphical  sketch  of  the  computed  values  is  given  in  Fig.  18. 

(Note  that  the  true  response  is  not  shown  in  Fig.  18). 


IV.  BILINEAR  AND  VOLTERRA  RESPONSES  OF 


1-M  SECTION  TO  SQUARE  PULSE  INPUTS 


Our  interest  in  the  response  of  nonlinear  Volterra  systems  to  square- 
pulses  arises  because  of  the  many  advantages  this  type  of  input  offers  in 
the  laboratory.  These  are:  (a)  it  is  easily  generated  in  the  laboratory, 
(b)  its  spectral  width  (f^)  is  readily  specifiable  through  the  pulse  width 
(b  =  1/f^),  and  (c)  the  formulas  for  the  system  response  are  tractable. 
Also,  as  we  shall  see  in  the  next  section,  it  is  possible  to  perform  para¬ 
meter  identification  via  pencil-of-function  method  for  such  input-output 
pairs. 

In  addition  to  the  Volterra  response  of  the  1-M  section  (see  Fig.  19), 
wc:  will  also  be  interested  in  its  bilinear  response,  because  the  bilinear 
response  V2{xj^,X2}  can  be  determined  in  the  laboratory  by  performing  tests 
upon  the  system  with  the  inputs  x^^  +  and  x^  -  and  then  using  the 
formula  ^ 


Alternatively,  one  can  test  the  system  with  the  inputs  x^^. 
Then  the  bilinear  response  is  given  by 


(37a) 

x^  and  x^  +  x^ . 


^2  2  [72^^!  ^2^  ~  ^^[x^]  -  V2[X2]J 


(37b) 


Fig.  19  General  1-M  Section 

2Formulas  (37a)  and  (37b)  assume  that  the  Kernel  112(3^^,82)  is  symmetric. 


Of  course,  the  re^l  reason  for  the  Interest  in  the  bilinear  response  Is 
Its  bilinearity  property.  For  example,  the  bilinear  response  to  two  square 
pulses  Pj^(t)  »  A(u(t)-u(t-Tj^))  =  A(u^-Uj^)  and  P2(t)  =  B(u(t-T^)-u(t-T)  )=B(Uj^-U2) 
Is  given  as 

“  AB[V2{u^,u^}  -  V2{uj^,u^}  -  V2{u^,U2}  +  V2{u^,U2}]  (38) 


4.1  Bilinear  Response  to  Step  Inputs 

4.1.1  Elementary  1-M  Section 

Consider  that  Xj^(t)  =  u(t  -  a)  and  X2(t)  =  u(t  -  B)  with  3  >  ct.  Then 
the  associated  bilinear  response  of  the  elementary  1-M  section  (see  Fig.  3) 
Is  given  by 


-asj^  -Bs^ 
e  e 


^2^®1’®2^  =  H2(Sj^,S2)  X2(s2) 

^  _ 1 _  _  _ 

(s^+a)  (s2+b)  (Sj^+s2+c)  s^^  S2 

11111  "^^2 

B  — i.(  ^ - i - )  (  1 - i - )  ^ ® - 

ab  Sj^  Sj^+a  S2  S2+b  s^  +  $2  +  c 

so  that,  by  application  of  George's  theorem  [8],  the  bilinear  response  becomes 


V2{Xj^(s),X2(s)}= 


ab  ^s(s+c) 


£!1  [-1  _  _ 

ab  ^cs  (c-a)  s+a 


(s+b) (s+c) 
1 


_  ^-a(B-a)^ 


(s+a) (s+c)  (s+a+b) (s+c) 


) 


(  c-b  )  s+b 


g-a(B-a) 

(c-a-b) 


s+a+b 


+  e 


-a(B-a) 


c-a 


1 

c-a-b 


ab  Lcs 


-a(B-a)  ^  -a(B-a)  . 

e _  .1  1 _  e _  1 

(c-a)  s+a  (c-b) (s+b)  (c-a-b)  s+a+b 


+  (- 
c 


_b _ 

(c-b) 


be-a(B-”)  1  -1 

(c-a) (c-a-b)  s+c  J 


(39) 
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Therefore  the  bilinear  response  in  time  domain  is 


V2{xi,X2}' 

where  L  = 

■  ;^[d^-d^Le"®‘''-d2e‘''’'’+d3Le"^®^'’^*''  +  (d^-d5L)e"‘'‘'']u(t') 

e-a(e-“).  t-  =  t  -  8  and 

(40a) 

= 

1/c 

(40b) 

'^l  = 

1/ (c-a) 

(40c) 

^2  = 

l/(c-b) 

(40d) 

= 

1/  (c-a-b) 

(40e) 

^  = 

b/c  (c-b) 

(40f) 

^  = 

b/  (c-a) (c-a-b) 

(40g) 

Similarly,  for  the  case  a  >  3,  we  obtain 

V{x.  ,x„}  =  [d  -d,e"^*^’-d.Me"*’’^'+d,Me“^^'^^^’^'  +  (d,-Md,)e"‘^'^' ]u(t')  (40h) 

izabol  L  j  0/ 


where  M  =  t'  =  t-a  and 

d^  =  a/c(c-a) 

d^  =  a/ (c-b) (c-a-b) 

Note  that  d,  -  d.  =  d,  -  d-. 

4  5  6  7 

4.1.2  General  1-M  Section 


(40i) 

(40j) 


For  the  general  1-M  section  of  Fig.  19,  the  bilinear  response  to  the  unit 
steps  Xj^(t)  =  u(t-a)  and  X2(t)  =  u(t-6)  is  obtained  immediately  from  formula 
(17)  and  (40a).  That  is, 

V{xj^,X2}  = 


.  .  ,  ,  a  h .  o  1  2 

i,j,k=l  i  j 


n 

Z 


c,  t' 


+  d^-^L-'e  4  5 
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where  the  constants  are  defined  as  in  (40)  with  a^,  and  substituted 
In  place  of  a,  b,  and  c.  For  example. 


jijk  . 


"k  -  ®1 


and 


Ijk 


-a^(&-a) 


4.1.3  Symmetric  1-M  Section 

Since  the  bilinear  response  can  be  measured  only  for  the  symmetric  Kernel 
It  is  useful  to  give  the  counterparts  of  (40a)  and  (41)  for  this  case.  For 
the  elementary  1-M  section,  setting  b  =  a  we  have 


V{xj^,x2}=  -j  [  d^-(d^L+d2)e“®’^’-Nl3Le~^^’^’  +  (d^-d^L)e”‘^’^']u(t') 


-ct' 


(42) 


Similarly,  for  the  general  1-M  section,  setting  b^  =  a^ ,  we  have 
V{xj^,X2}  = 

?  _J_  [  .  ,ljk  ^ljk^-^l^’.,ljk 

i.j.k=.i  Vj  °  ^  2 


+  di^j^  e  ^  i  +  (d^J^  -  d^J'^  L^^^)e  ]  u(t') 

(43) 


-(a^+a3)t'  -c,.t’ 

*3  “  '2  ■  '-4  -5 

4.2  Bilinear  Response  to  Square  Pulses 
4.2.1  Elementary  1-M  Section 

We  are  interested  in  finding  the  bilinear  response  of  the  1-M  section 
to  the  input  pair  p^(t)  =  u(t)  -  u(t-Tj^) ,  P2(t)  =  u(t-Tj^)  -  u(t-T).  The 
associated  bilinear  response  (in  s-domain)  is 


^(2)^®1’®2^  =  ^2^®1’®2^  X(S3)X(S2) 


(Sj^+a)  (s2+b)  (S3+S2+C)  Sj^ 


-T  s  -T  s  -Ts, 
(1-e  ^  ^)  (e  ^  -e  ^) 
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•1 


Fig.  20  Square  pulse  -  Input  pair 


From  this  we  can  obtain  the  bilinear  response  by  application  of  George's 
theorem.  Alternatively,  the  desired  bilinear  response  is  computed  using  the 
bilinearity  property: 


^2W’P2^  =  V2tu^,u^}  -  V2{uj.u^}  -  V2{u^.U2}  + 

where  u  =  u(t),  u,  =  u(t-T, )  and  u„  =  u(t-T). 
o  112 

By  use  of  formula  (AO) ,  we  obtain 

'^2^Pl’P2^  = 

—  [d  -d  d-e~'^‘^'+  d,L,e"^^'^^^*^’  +  (d,-d,.L,)e"‘^'^' ]  u(t’) 

ab  oil  2  31  451 


r.  ,  -at'  ,  -bt'.j  -(a‘H))t'..,  ,  .  -ct' .  , 

-  [d  -d,e  -d„e  +d-e  +{d.-d,.)e  ]u(t  ) 

o  1  2  3  4  5 


-  [d  -d  L  d-e  +  d  L  e  +  (d  -d  L  ]  u(t") 

o  1  2  3  4  5 

+  [d  -d  L  d-e'^'^'V  d,L  (d  -d  ]  u(t") 

ol2  2  32  452 


(44a) 


t'  =  t  - 
t"  =  t  -  T 


(44b) 

(44c) 

(44d) 

(44e) 

(44f) 
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Let  us  also  define,  for  convenience 


-bT. 


■cT„ 


Then  simplification  of  (44a)  yields 


f  0 
(1-Lj 


1^  -at*  j  -(a+b)t'  .  .  -ct',  /. ,s 

— [d^e  -d^  e  '  +  d^  e  ]  u(t') 


I  ^  Cdj^d-Mj)  =-<‘*«''Vd5(N2-Lpe-“"l 


(44g) 

(44h) 


for  0  _<  t  < 
for  £  t  <  T 

for  T  £  t 

(441) 


4.2.2  General  1-M  Section 

For  the  general  1-M  section  of  Fig.  19,  the  bilinear  response  to  the 
square  pulses  Pj^(t)  and  P2(t)  Is  obtained  Immediately  from  formulas  (17) 
and  (441).  That  Is, 


f  0 


for  0  <  t  <  T, 


(1-L,^J^) 


41  a,  t  ,  , .  “  (a ,  +b . )  t  ,  , ,  *"*“1,  ^ 

l  E  ABC  - i e  ^  -d^^^  e  ^  ^  +d^^*‘  e  ]  u(t') 

(  l,J,k-l  ‘  J  ''  •l’’!  ‘  5  5 

'  for  <  t  <  T 


l,j,k-l  ^  J  Vj  3  2  2  5  2  2 


(a,+b,)t" 


V 


for  T  <  t 


where  the  constants  are  defined  as  In  (40)  and  (44)  with  a. ,b.  and  c, 

i  j  K 

substituted  In  place  of  a,b,c.  For  example 
dJJ*^  -  l/(c^-a^) 


and 


,ijk  -Vi 

^1  “  ® 


lu(t") 

(45) 
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4.2.3  Symnetric  1-M  Case 


Since  the  bilinear  response  can  be  measured  only  for  the  symmetric  Kernel, 
it  is  useful  to  give  the  counterparts  of  (441)  and  (45)  for  this  case.  For 
the  elementary  1-M  section,  setting  b  =  a  we  have 

fo  for  0  <  t  < 


^  ^1  r  .  -at'  ,  -2at'  ,  ,  -ct', 

/  — 2~  I  ^  -  d^  e  +  d^  e  ]  u(t') 

\  a 


for  <  t  <  T 


(46) 


1-L 


^  [  d3L2(l-L2)e"^^*^'  +  (N2-L2)d'‘^*'"]  u(t") 

\  a 

for  T  <  t 

Similarly  for  the  general  1-M  section,  setting  b^  =  a^ ,  we  have 


V**1’P2^  “ 


(  0 


1-L 


ijk 


< 


A,B.C. 


k-l  ‘  J  >' 


for  0  <  t  <  T, 


-a,t'  -(a,+a.)t'  -c,  t' 

1  -dUk,  ‘  3  +dl3ke  k  ]  u(t') 

for  <  t  <  T 


n 

E 


(l-L^J*^) 


"|  u(,") 


A.B.C, 


3  2 


Consider  an  elementary  1-M  section  with 


H^(s)  »  H^^(s)  =  ^  ^ 


H^(s)  =  4/(s+12) 

and  p^(t)  =  u(t)  -  u(t-l),  p2(t)  =  u(t-l)  -  u(t-2) 


Then  using  (4)  we  have 

V{pj^,P2}  = 

0 


( 


•l*’! 

121 

__ 


for  0  <  t  <  T 


111  ^1*^'  111  (ai+b^)t'  2^  -c  t' 

d^  e  -d^  e  +  d^  e 


121  -aj^t’  2J.  -(a^+b2)t'  ^21 

dj^  e  -dj  e  +  d^  e 


211 

^2^ 


221 

A2B2Cj^(1-L^  •^) 
^2^2 


",  -a„t'  211  "^®2'^1^‘^’  211  ■'^l’^' 

,211  2  e  +  e 

d^  e  3  5 


221  ^2^'  221  (a2'^2^’^'  ,221  “'^1*^' 

dj^  e  -d^  e  +  d^  e 


) 


u(t') 


for  £  t  <  T 


,111  .111  ..  Ill-  ■^^I'^^l^’^",  .„111  111,  ■'"l*'' 
^3  ^2  ^  ®  +  (N2  -L2  )e 


121  121  121-  ^121  121,  "'^1*'' 


211 

A^B^C^d-L^  n 
^?i 


211  211  211-  . ,211  -211,  ~‘^l‘^' 

“3  ^2  "^^^2  ~^2 
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221  I— 
A2B2Cj^(1-L^  ) 

- - -  jj 

^2*^2  L 


221,221,,  ,.221,  ",  ,„221  ,  221,  “‘^l 


d3“L2“(l-M2  )e 


where 


h  =  1 


=  Bi  =  1 


Cl  =  12 


Cl  =  A 


T  =  2 


+(N2  -L2  )e 


for  T  <  t 


C2  =  0 


S  =  0 


u(t") 


V{pj^,P2)  = 


-1,  1  -t'  1  -2t'  1  -12t' 

4(l-e  )^-e  --e  +—  e 

.  -40(l-e“^)  r  1  -t'  1  -6t'.  5  -12t'l 

(5)  [ir  ®  “  6  ®  J 

,  -40(l-e'^)  1  -5t'  1  -6t’  .  1  -12t'l 

^  - Tsl -  L  7  "  -  6  ®  "  J 


,  -40(l-e  ■")  1  -5t'  1  -6t 

^  Tsl  7"  "6® 


400(l-e 


'1  [i 

7 


-5t'  1  -lOt'  5  -12t' 

2  -  TT  e  +  ^  e 


for  0  <  t  <  1 


u(t’) 


for  1  <  t  <  2 


{L(^  -h  Ll  -In  r  -12  -U  -12t"  ' 

^vl-e  )  Yo  ®  (1-e  ;e  +(e  -e  )e 

-40(l-e~^)  1-1,,  -5x  -6t". ,  -12  -1,  -12t" 

+  (l-e  )e  +(e  -e  )e 

-40(l-e~^)  fl  -5.,  -1^  -6t''  ,  -12  -5,  -12t"l 

+  (5)  0  e  (l-e  )e  +(e  -e  )e 

400(l-e~^)  [  1  -5,  -5,  -lot"  .  -12  -5,  -12t"l  \ 

^  ^25)  [2  ®  J  / 


-40(1-£ 


u(t") 


for  2  <  t 
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0 


for  0  <  t  <  1 


^0,2298626"^*^"^^  -  0.252848e“^^*^  0. 022986e"^^ 

-0.459724e~^‘^"^^  +  0.8428276”^^*^"^^  -  0.3831036"^^^*^"^' 

-1.135157e"^^*^“^^  +  1.324349e"^^*^"^^  -  0. 189193e"^^^‘^"’-^ 


2.270313e  -  7.9460960"^°^*^  +  5.6758  e  !•  u(t-l) 


( 


for  1  <  t  <  2 


) 


0. 0587980  -  0.9301610“^^^^"^^ 


-0.3079690"^^*^"^^  +  1.8603220"^^^’^"^^ 


-0.0056410 


-6(t-2)  0.0534920”^^^’^"^^ 


+0.0531790 


-10{t-2)  -  0.1069830 


-12(t-2) 


} 


u(t-2) 


for  2  <  t 


{ 


for  0  <  t  <  1 
-2(t-l) 


-O.2298620"^*^~^^  +  1.1351560"^^^  -  0.2528480 


+  2.1671760“^^"“^^-  7.9460960"^°^*'  ^^+  5.1265  0 


->}  „ 


(t- 


for  1  <  t  <  2 


^  0.0587980  O,3136l0~^^*^"^^ 

+  0.0531790'^°^'^“^^+  O.876670"^^^’^"^^  |  u(t-2) 

for  2  <  t 
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4. 3  Volterra  Respoo&e  to  Step  Input 
4.3.1  Elementary  1-M  Section 

The  Volterra  response  to  a  unit  step  u(t)  is  obtained  by  setting 
a  “  S  =  0  and  L  =  1  in  (40a) : 


y2(t)  =  ^'^i-d^e  ^^'*’*’^*^+(d^-d^)e  ‘^‘^]u(t)  (48) 

where  d  .  d^  are  the  same  as  in  (40). 

o  3 

4.3.2  General  1-M  Section 

For  the  general  1-M  section  of  Fig.  19  the  Volterra  response  to  a  unit 

i  "i  Ic 

step  u(t)  is  obtained  by  setting  a  =  g  =  0  and  L  =  1  in  (41): 
y2(t)  = 


n 

Z 


1  r  .iik  .lik  ^i*"  ,iik  -b^  t 
— —  [d  J  -  d/  e  -d  e  J 
a ,b .  o  1  2 


i,j ,k=l  i  j 

+  d^j*^  (dJJ*^-  u(t) 

4.3.3  Symmetric  1-M  Section 

For  the  symmetric  elementary  1-M  section  we  have  from  (48) 


(49) 


y2(t)=  ^ld^-(d^+d2)e~^*'+  d^e'^^V  (d^-d^)e"^’']  u(t) 


(50 


Similarly,  for  the  multiple  pole  case  we  have  from  (49) 


y2(t)  = 


n 

Z 

l.J,k=l 


1 

l"i 


.  ,ljk  jljk  1  jijk  J 
[  d  -  d,  e  -  d„  e 
o  1  Z 


.  -(a,+a.)t  ...  ,  -c,  t 

+  d^J*^  e  ^  1  +  (dJJ*^-  d^^Se 


4 . 4  Volterra  Response  to  Square  Pulse 
4.4.1  Elementary  1-M  Section 

The  Volterra  response  to  a  square  pulse  p(t)  = 


u(t) 


u(t)  -  u(t-T)  A  u  -  u,  is 
—  o  i 


(51) 
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By  use  of  (AOa)  and  (40h),  and  some  simplification,  we  have 

y2(t)  - 


d,  d 


-bt 


+  d. 


e-(a+b)t^ 


(d, 


-ct'. 


]  u(t') 


0  <  t  <  T 


;^[d3(L-l)(M-l)  e  ^^^^*'+((d^-d5)N+d5L+d^M-d5-d^)e"‘'’']  u(t) 

T  <  t  (52) 

^  c^r 

where  t'  ■  t-T,  L  •  e  ,  M  =  e  and  N  =  e  .  All  other  symbols  are  the 
same  as  in  (40) . 


4.4.2  General  1-M  Section 

For  the  general  1-M  section  of  Fig.  19,  the  Volterra  response  to  a 
unit  pulse  p(t)  ■  u(t)  -  u(t-T)  is  obtained  immediately  from  formulas  (17) 
and  (52) .  That  is 
y2(t)  - 


?  Vii  ijk  _  ^ijk ,  .  ,ijk 

i.j.k-i  “I'-j  ^  2 

.  ,ljk  ,,ljk  ,ijk.  . 

+  d^-^  e  -^  +  (d^-J  -  d^-^  )  e  ]  u(t) 

for  0  <  t  <  T 


-c,  t' 


((d^J^-  d  d  d^^*^)e  ^  ]  u(e) 

4  5  5  7  5  6 


for  T  <  t  (53) 

4.4.3  Symmetric  1-M  Section 

The  counterparts  of  (52)  and  (53)  for  the  symmetric  case  can  be  written 
immediately,  but  are  not  given  for  brevity. 
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Example  10 


H  (s) 
a 


"  s+12 


Then  using  eq.  53  we  have 
V{pj^,P2) 

r 

1 


{ 


,  .  ^  1  -t  I  -t  1  -2t^,  1  _ I.  -12t, 

12  "  TI^  n®  10®  "^^132  110^®  ^ 


40,  1  1  -t  1  -5t_^  1  -6t_^,  5  5,  -12t, 

-  12  -  ir  -  7®  ^  -  66>®  ^ 


40,  1  1  -t  _  J_  g-t  +  1  g-6t  ,_1 _ 1_.  -I2t 

-^[  Y2  ■  7®  11  6  U32  42-’®  J 


11  ''  '6 

+  i  1  -5t.  1  -lot 

25^2  “  7 


,  ,  5  5.  -12t, 

■"^84  ••  14>®  1 


! 


u(t) 


for  0  <  t  <  1 


{,,  1,-1  1,2  -2t’  ,,  1  1,  -12,  1  -1. _ 1-1  _ 1 

4[io(e  -1)  e  +(<132  "  no^®  II^  TiO^  110 


li  -12t', 

132^®  ^ 


-6t' 


3s  -12_,  5  1-5 


-  ^[|(e-'-l)(e-5-l)e—  -l((3^  -  -®  42' 


66  132''®  ' 


4U,i,  -i  1,,  -5  1,  -6t'  ,,  1  Is  -12  1-5 

“  “5^®  ■^^®  ■^^^r32  ■  «^®  ■^42® 


40,1,  ^-1_ 

5-11  5s  -12t'  ] 


r)e 


66  42  84' 

*  if  - 

for  1  <  t 
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[0.3333  +  0.7273e~*^  +  0.4000e~^’^-  2.2856”^*^ 
-2.6667e"^*^  +  S.OOOe”^®*^  -  4.504e"^^*^]  u(t) 

for  0  <  t  <  1 

[0.1598e"^^’^~^^  -  1.6743e~^^*^"^^  +  7. 8926e“^°^*^"^^ 
-5.7455e"^^^‘^‘^^]  u(t-l) 

for  1  <  t 
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V.  IDENTIFICATION  OF  H^Cs^.s^)  USING  PULSE  INPUTS 

In  this  section  we  shall  deal  with  the  central  problem  of  the  report, 
i.e.,  identification  of  H2(Sj^,S2).  It  will  be  shown  that  this  can  be  accom- 
plished  by  use  of  either  the  bilinear  response  to  a  pulse  pair  ^^(t)  = 
u(t)  -  u(t-Tj^)  and  P2(t)  *  u(t-Tj)  -  u(t-Tj^-T2),  or  the  Volterra  response  to 
a  square  pulse  p(t)  =  u(t)  -  u(t-T)-  It  is  assumed  that  the  transfer  function 
H2(s^,S2)  is  symmetric  or,  equivalently,  that  H^(s)  =  H^(s)  in  the  block  dia¬ 
gram  of  (19).  This  assumption  is  entirely  harmless  because  one  can  only  identify 
the  symmetric  equivalent  a  Volterra  system  from  input-output  measurements  [  7  ]. 

The  problem  is  considered  in  two  parts.  First  the  estimation  of  poles 
of  and  is  considered.  Next,  the  residues  are  estimated. 

5 . 1  Identification  of  Poles 

5.1.1  Poles  from  Bilinear  Response 

For  a  symmetric  1-M  section  the  bilinear  response  to  the  square-pulse 

pair  p^(t)  =  u(t)  -  u(t-T^)  and  ~  u(t)  -  u(t-T)  was  shown  to  be  given 

by  (47).  Over  the  time  interval  T^^  <  t  <  T  this  response  can  be  visualized  as  the 

step  response  of  an  equivalent  linear  system  shown  in  Fig.  21,  where  a^jAa^+a^  and 

the  residues  P  ,  Q  and  R  defined  according  to  the  second  part  of  (47). 

i  ij  1 

Clearly,  the  pencil-of-functions  method( 10] , ( 12] (.3^  be  used  to  find 
a^  i  =  1,  . . . . ,  n 

,i  =  1,  . .  ,i 

c^^  k  =  1 . m 

Note  that  in  general  the  total  number  of  poles  is 


N  =  n  +  + 


(54) 


If  the  1-M  section  is  completely  symmetric,  i.e.,  if  =  H^, 

then  the  total  number  of  poles  is  again  given  by  (54);  however  the  poles 


a^  occur  with  a  multiplicity  of  2,  i.e.. 


^The  reason  for  the  Interest  in  the  bilinear  response  is  the  extra  control  on 
the  energies  of  the  various  modes  provided  by  the  widths  of  the  pulses  Pj^  and 

P2- 
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{T^  for  the  response  over  £  t  <  T 

T  for  the  response  over  T  £  t 

Fig.  21  Equivalent  linear  system  for  bilinear  response 
of 

We  remark  that  the  type  of  H2(sj^,S2)  considered  by  Ewen  were  completely 
symmetric  1-M  type.  The  simulation  examples  given  below  all  pertain  to  com¬ 
pletely  symmetric  1-M  sections. 

Example  11 

Consider  a  1-M  section  with 


50.5 


H^(s)  -  H.  (s)  -  H  (s)  =  - 

s  +S+25. 25 


The  bilinear  response  to  square  pulses  Pj»P2  with  T^^  =  1  s.  and  T2  =  1  s. 

(so  that  T  “  2  8.)  was  generated  over  the  interval  T^  £  t  <  T. 

Using  a  sampling  Interval  A  =  0.02  s.,  the  response  was  identified 
using  the  computer  program  IGRAM.  The  following  s-domain  poles  were  obtained 
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-0.553  -  j5.023 


=  -0.553  +  j5.023 

=  -1.135 +  j4. 9792 
=  -0.998  -  J10.025 

=  -0.998  +  jlO.025 

Sg  =  -1.1396  -  j9.9968 

=  -1.0021 

The  normalized  mean  square  error  was  v  ■■  -0.0012.  Note  that  the  following 
program  parameters  were  used  IREM  =  0,  IDLY  =  0. 

By  observation  it  is  deduced  that 
a^^  =  0.553  +  j5.023 

a.^  =  0.553  -  j5.023 


Over  the  time  interval  T  j<  t  the  bilinear  response  is  given  by  the 
third  part  of  (47),  and  can  be  visualized  as  the  unit  impulse  response  of 
a  linear  system  as  shown  in  Fig.  22. 

Again  the  pencil-of-functlons  method  can  be  applied  to  determine 
a^^  i=l,  ....,n  j“l,  ....,  i 

. ® 


Note  that  the  total  number  of  poles  is 


n(n+l) 
N  =  — 2 - 


(55) 
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y,(t-T) 


Fig.  22.  Equivalent  linear  system  for  quadratic  response  of  H2(sj,s^) 


Note  that  a..  =  a.  +  a..  Valid  for  T  <  t 
U  I  J 


5.1.2  Poles  from  Volterra  Response 

For  a  1-M  section  the  Volterra  response  to  the  square  pulse  p(t)  = 
u(t)  -  u(t-T)  was  shown  to  be  given  by  the  first  part  of  (53).  For  the 
symmetric  case  this  response  can  be  shown  to  equal  the  step  response  of 
an  equivalent  linear  system.  This  linear  system  has  the  same  form  as  in 
Fig.  21  with  T  =  0.  Recall  that  a^^  =  a^  +  a^  and  the  residues  are  defined 
suitably. 

Again,  the  pencil-of-f unctions  method  can  be  applied  to  determine 


1  Tl 

j  =  1.  ,  1 

‘Tc 

Ic  "  ly  in 

In  general,  the  total  number  of  poles  is 


.  n(n+l)  , 
N  =  n  +  — 2 - 


(56) 


If  the  quadratic  is  completely  symmetric,  i.e.,  if  =  H^,  then  the 
total  number  of  poles  is  again  given  by  (56);  however,  the  poles  a^  occur 
with  a  multiplicity  of  2. 
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Over  the  time  Interval  T  ^  t  the  quadratic  response  is  given  by  the 
second  part  of  (53),  and  can  be  visualized  as  the  unit  Impulse  response  of 
a  linear  system  as  shown  in  Fig.  22. 

Again  the  pencil-of-functions  method  can  be  applied  to  determine 

®ij  i  =  1 . .  j  =  1 . i 

Cj^  k  =  1, . . . .  ,m 


Note  that  the  total  number  of  poles  is 


^  n(n+l) 
2 


(57) 


5.1.3  Remarks 

1.  The  test  engineer  has  a  choice  here  between  the  use  of  the  two  segments 

of  Cor  y2(t))«  At  this  time  it  appears  that  the  use  of  last  segment  is 

preferrable,  since  the  dimensionality  of  identification  is  lower  in  this  case 

(as  evidenced  from  a  comparison  of  (54)  with  (55)),  although  quite  extensive 
experimentation  is  necessary  to  make  a  definitive  statement. 

2.  Since  the  identified  values  a. .  and  c,  will  not  coincide  with  the 

ij  k 

true  values,  it  is  necessary  to  isolate  the  poles  by  visual  inspection  or  by 
a  suitable  computer  routine.  For  example,  if  these  numbers  for  the  case 
(n,m)  =  (2,1)  are 

2.1,  4.15,  9.8,  5.95 


then 


*11 


12 


22 


2.1 


4.15 

5.95 


9.8 
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and  we  may  take 


a^^  =  1.05 


3.1 


9.8 


5.2  Identification  of  Residues 


After  the  poles  have  been  determined,  we  can  write 


i 

E 


m 


Vj^k 


(58a) 


where  are  defined  in  accordance  with  (47)  or  (53)  (these  are  known 

functions  of  the  poles  and  the  pulse  width  T) .  Further,  by  defining 


e  =  A.A.C, 
V  1  J  k 


V  =  i  +  j+  k-  2 
y  =  0, . . . ,M-1 


(58b) 

(58c) 


we  obtain  the  following  set  of  simultaneous  equations 


Y  -  G  E  (59) 

where 

Y  =  [y2(0)  y2(^) . y2((M-l)A)]^ 

G  “  [8y(p)]  M  X  I  matrix 

E  =  [e^  .. .  e^] 


Note  that  M  denotes  the  number  of  time  samples  used  in  (58)  and  (59), 

and  I  denotes  the  nunfcer  of  unknown  residue-products.  If  M  is  chosen  equal 


If  some  of  the  poles  are  complex  (occurring,  in  conjugate  pairs),  the 
associated  residues  are  also  complex.  In  such  cases,  it  is  possible 
to  equate  real  (and  imaginary)  parts  on  both  sides  of  the  equation  (59)  to 
obtain  real  coefficient  equations 


53 


to  I  then,  we  obtain  the  solution 


G“^  Y 


(60) 


More  generally  If  M  >  I,  then  we  obtain  the  least-squares  solution  [15] 

(61) 

Finally,  the  equations  (58b)  can  be  solved  straightforwardly  to  yield  the 


E  =  (G’^G)"^  G'^y 


residues  and  Cj^.  Because  of  the  homogeneity  of  the  relations  any  one  of 


the  real  residues  (or  the  magnitude  of  a  complex  residue)  can  be  taken  as  1. 


VI.  EXAMPLES 


Two  computer-generated  examples  will  be  presented.  The  first  is  a 
simple  case  with  n=m=l  and  is  intended  to  clearly  present  the  details  of 
the  procedure.  The  second  is  a  more  complicated  case,  representing  a  some¬ 
what  realistic  channel.  It  has  a  linear  transfer  function  characterized 
by  a  6th  order  butterworth  filter,  and  a  quadratic  transfer  function  with 
(m,n)=(3,l) . 

Example  12 

Consider  that  the  response  of  a  quadratic  subsystem,  with  H  (s)  = 

l/(s+l)  and  Hj^(s)=l/ (s+0. 5) ,  to  a  square  pulse  of  one  second  duration  has 
been  measured.  With  a  sampling  interval  A=  0.02,  q  =  0.8,  and  input  a 
mit  sample  pulse  (digital  impulse) ,  the  symmetric  Gram  matrix  of  the  in¬ 
formation  signals  turns  out  to  be  (only  lower  triangular  entries  are  shown) 
~0.88303D+01  ~ 

0.410580+02  0.195070+03 

0.188430+03  0.908280+03  0.428710+04 
-0.16677I1+01  -0.463240+01  -0. 128670+02  0.277780+01 

_-0. 913830+01  -0.336160+02  -0.116230+03  0.771600+01  0.351510+02_ 

The  normalized  cofactor  square-roots  (y^^//D^)  are 

1.0000  0.35933  0.032129  0.22357  0.026502 

5 

By  use  of  (6)  and  (7b)  the  negatives  of  the  poles  are  computed  to  be 


^  The  reason  to  associate  1.9997  with  a^^^^  is  that  when  segment  1  is  analyzed, 
pole  at  0.9996  also  shows  up. 
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1.9997 


^11 

=•  0.5002 

so  that 

=  0.9998 
=  0.5002 

Note  that  the  waveform-fit  error  in  modeling  the  poles  was  almost  negligible 
(normalized  mean-square  error  =  10  pointing  to  the  success  of  pulse  testing 
approach. 

Now,  using  these  poles  we  find 

y2(t')  =  (-0.2664375  '+0. 4119195e"°‘ ) 

f  6 

Using  a  single  point,  y2(t  -0)  =  0.292924,  we  obtain  =  2.013,  so  that 

A^=  1 
^1  =  2.013 


Example  13 

Linear  TF  Hj^(s)  Sixth  order  butterworth  with  cutoff  f^  =  10  MHz 


Hj^(s) 


[A' 


6.1528908(10^^) 

+  2.4276(10®)s^  +  2.9467(10^^)8^  +  2.2676(10^^)s^  + 

T9  9  TQ  46 

1.1633(10  )s  +  3.78358(10  )s  +  6.1528908(10  )] 

6.1528908(10^^) _ 

+  242.761^  +  29467X‘‘  +  2267581X^  +  1.6133(10®)X^ 

q  10 

+  3.78358(10  )A  +  6.1528908(10  )] 


A  better  value  A^A^^Cj^ 


2.0007  is  obtained  using  50  points  and  formula  (21). 
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Where  A  =  (10  )s  (units  of  Mrad./s) 


50.5(10^) 


Quadratic  TF  H  (s)=H  (s)=  — r - = - 

^  +  (lO^s  +  25.25  (lO^S 


X  +  lOX  +  2525 

The  above  2-variable  system  was  excited  by  a  square  pulse  p(t)  of 
duration  T  “  1  sec  and  the  response  y^(t)  recorded.  The  system  was  next  excited 
by  the  opposite  polarity  pulse  -p(t)  and  the  corresponding  response  y  (t) 
also  recorded.  From  these  responses  we  obtained  the  linear  response  yj^(t) 
and  the  quadratic  response  y2(t)  by  use  of  (11). 

The  results  of  identification  are  given  below: 

Estimated  Linear  Transfer  Function  - 


15.2  X  +  6.1523(10  ) 

tX^  +  242. 76X^  +  29467X^  +  2267593X^  +  1. 6133(10®)X^  +  3.784(10^)X 


+  6.15237(10^  )  ) 


Estimated  Quadratic  Transfer  Function 


^  ^  0.003  X  +  505.07 

H  (X)=  H  (X)=  - 

^  X  +9.9998X  +  2524.8 


^  It  is  easy  to  show  that  yj^(t)  =  ■|^[y  (t)  -  y  (t)],  and  y2(t)=  (t)+y  (t) 
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APPENDIX  A 


GEORGE'S  THEOREM 


This  theorem  helps  evaluation  of  the  quadratic  volterra  response  (or 
the  bilinear  response)  by  inspection  from  the  associated  response. 

Theorem  If 


Y(2)(Si,S2) 


(Sj^+a)  (s2+b) 


C(s^+ 


S2) 


(Bl) 


then 


Y2(s)  = 


C(s) 


s  +  a  +  b 
Proof;  It  is  readily  shown  that 

joo 


Y2(s) 


1 

2TTj 


C(s) 


(Sj^+a)  (s-s^+b) 


dsi 


(B2) 


(B3) 


(see  Schetzen  [7]).  Then  the  desired  result  follows  Immediately  by  use  of 
the  residue  theorem. 


Corollary  Let 

-asj^  -6S2 


^2  ^®1»^2^  “  (s^+a)(s2+b) 

(B4) 

If  6  >  a,  then 

Y  -  -a(6-a)  e  ^®C(s) 

2  (s  +  a  +  b) 

(B5) 

If  a  >  B  then 

Y  (s)  =  e-<""c(s) 

2''®'’  ®  (s  +  a  +  b) 

(B6) 

